Calculation of transport coefficient profiles in modulation experiments as an inverse 

problem 
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The calculation of transport profiles from experimental measurements belongs in the category of 
inverse problems which are known to come with issues of ill-conditioning or singularity. A reformula- 
tion of the calculation, the matricial approach, is proposed for periodically modulated experiments, 
within the context of the standard advection-diffusion model where these issues are related to the 
vanishing of the determinant of a 2x2 matrix. This sheds light on the accuracy of calculations with 
transport codes, and provides a path for a more precise assessment of the profiles and of the related 
uncertainty. 
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Inverse problems, i. e. going from data to model 
parameters, are ubiquitous in science and engineering. 
They are known to come with issues of ill-conditioning 
or singularity, which make difficult relating the accuracy 
of the computed model parameters to the errors in the 
input data. Inferring transport properties of magneti- 
cally confined plasmas from perturbative experiments [l[ 
belongs in this category. This issue is crucial to both the 
theoretical understanding of transport processes and the 
practical control of fusion plasmas. Classical reconstruc- 
tions of transport coefficients rely on standard transport 
codes, which ignore the possible ill-posedness of the prob- 
lem by employing ad hoc regularizations. To the best of 
our knowledge, so far the only authors that took it into 
account were Andreev and Kasyanova [H, who provided 
a detailed analysis of the uncertainties present in the 
case of a localized impulsive source. However, as shown 
later, experimentalists happened unwittingly to avoid ill- 
conditioning and singularity by trying naturally to sepa- 
rate the domain where transport is measured from that 
where sources are present. Unfortunately this separation 
is almost impossible for measuring the transport of angu- 
lar momentum, and only partially possible for heat trans- 
port. It is therefore most important to warn the users of 
transport codes, and to benefit from an alternative safer 
approach. For the first time this paper proposes such an 
alternative. The corresponding technique is elementary, 
quite general, and requires much less numerical calcula- 
tions than with transport codes. It may be applied to any 
experiment aiming at the calculation of transport coeffi- 
cients, also in media other than plasmas. Its only limita- 
tions are: (i) it works in effectively one-dimensional ge- 
ometries; (ii) applies to perturbative regimes, i.e., where 
relevant equations may be linearized around unperturbed 
equilibria; (iii) the forcing term causing the perturbation 
must be periodic in time. 

The standard procedure for inferring transport coeffi- 
cients is through solving the advection- diffusion for the 
generic quantity £(r, t) (which may stand for the pertur- 



bation to particle density, temperature, momentum,...): 
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In O Xi V are functions taken out of a set of trial profiles 
guessed a-priori, usually simple analytical expressions 
containing a certain number of free coefficients which are 
given explicit values by minimizing under some norm the 
gap between the experimental £ profiles and the recon- 
structed ones. The regularization built in this approach 
is straightforward: it consists in choosing well-behaved 
trial functions for Xi V. 

In this work we address the problem of solving Eq. ([1} as 
an inverse problem in the case of periodically modulated 
perturbations. We present a simple algebraic method to 
recover exact solutions; that is, x, V , may be exactly and 
directly (not by means of iterative procedures) obtained 
from the smoothed data without the recourse to any ad- 
justable trial function and minimization procedure. We 
identify under which conditions the problem is well- or 
ill-posed, i.e., the solution (x, V) does sensitively depend 
from the input data. We discuss how error bars in the raw 
measurements are propagated to final estimates of V). 
This is a delicate issue when Eq. (JTJ) is addressed as a 
least-squares problem using smooth trial functions for 
(x, V): from the one hand, they prevent any oscillation 
to blow up, thereby constraining the apparent error bars 
to values smaller than ought be. On the other hand the 
opposite limit is also possible: if trial functions are chosen 
that cannot match the true (x, V) profiles, the minimiza- 
tion procedure necessarily converges towards just a local 
minimum that is far away from the true global minimum. 
In this case, the error on the final (%, V) values is larger 
than expected from the raw data alone. Finally, we point 
out how data from different experiments can be combined 
to reduce the effective margin of error and thus yield a 
more accurate final estimate of (x, V). 
We consider a case with cylindrical symmetry and a 
purely sinusoidal forcing term. Calculations are easier 
using complex notation, thus the temporal behaviour is 
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factored out in the term exp(— itot) and Eq. (JTJ writes 

- iu( = r- x d r [r (xd r C - VQ] + S (2) 

After rearranging, integrating over the radial coordinates 
and using the boundary condition T(r = 0) = we get 



xd r ( ~V(= -r~ 



dz z (iuiC, + S) 



(3) 



The previous equation is reduced to an algebraic system 
of two real equations by expressing the signal in terms of 
amplitude and phase, ( = Ae 1 ^, and S = S r + iSf. 
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- J dz z (S r (z) — uiA(z) sin <j>(z)) 
o 

r 

i J dzz(Si(z)+ uA(z) cos 6(z)) 

(4) 

where the primes stand for differentiation with respect 
to r. Provided T(r) is known, Y(r) can be computed by 
inverting the matrix M(r). For this reason, in the follow- 
ing we refer to this method as the matricial approach. It 
should be mentioned that considering just time-periodic 
quantities in ([T]) leads to a considerable simplification: 
this reduces Eq. ([TJ from a partial differential equation 
to an ordinary differential equation ©. Several efforts 
were devoted to solve Eq. © using simpler techniques 
than its numerical integration, possibly even analytical 
methods 

We now solve Eq. (|4]) for Y(r) in terms of the eigenvalues 
A(r) and eigenvectors E(r) of the matrix M(r) 



Y = j/qEo + yiEi, yj 



A7 x r • E, 



U = 0, 1) (5) 



A solution exists wherever the matrix M is invertiblc, i. 
c. where none of the Aj's vanishes, which would imply 



det(M) = A 2 6' = 



(6) 



Excluding exceptional cases where A — 0, the singular 
points arc those where (/)' = 0. The origin is one such 
point, but there T = 0, too; thus ultimately r = is 
a regular point. Inspection of literature (see e.g. 
shows instead that the condition 6' = is usually ful- 
filled close to the location of the source. Qualitatively, 
this can be justified by noticing that, close to the source, 
the dynamics of C is ruled in Eq. more by the source 
than by transport if to is large enough, hence C « iSui -1 
and S' = implies 6' = 0. 

Physically, transport coefficients arc defined even at the 
singular points r = r s . Throughout this work we con- 
ventionally choose Ao to be the eigenvalue that vanishes 
at r s . Accordingly, the actual value of C produced by 



the source S must align the flux T exactly along the Ei 
eigenvector: 

M-Y = A 1 j/ 1 E 1 , r = 7l E 1 , (r = r s ) (7) 

Due to unavoidable uncertainties arising in the exper- 
iment, the estimated T takes on a component along Eq 
too, thus making the inversion impossible. However we 
may solve for the component of Y along Ei, whereas 
the component along Eo remains completely unknown: 
Y(r s ) belongs to a straight line parallel to Eo in the 
plane (x,V). 

We present now a comprehensive discussion, including 
both regular and singular points. We label with the sub- 
script "m" the quantities measured by the experiment: 
M m ,Cm,r m . Equation (0} implies M m • Y m = r,„. 
Let the subscript "c" label likewise the same quanti- 
ties as estimated by transport codes solving Eq. ([1]) 
via least squares. Finally, let "*" label the "true" (un- 
known) quantities that the experimental measure is ad- 
dressing: M* , C* , r* , and M* • Y* = IV T rn gener- 
ally contains some arbitrariness due to the lack of pre- 
cise knowledge about the source term, too. Let 5( = 
C c - Cm , SM = M C - M rn , 5*Y = Y c — Y* . The 
target of any transport modeling is 5»Y = 0. Finally, 
we write T c = T» + Ar s + Ar m + Ar c . In this ex- 
pression ATs is the error in the calculation of the flux 
due to the imperfect knowledge of the source, Ar m is 
the error related to imprecise measurement of £: £ m — £*, 
and the third term comes from the error in the recon- 
struction of the measured C : F° r brevity, we set 

r A = r* + Ar s + Ar TO , r c = r A + Ar c . Since Q is 

a first integral of (QJ, M c • Y c = T c holds. Using above 
definitions, we rewrite the latter expression as 

M m -tf*Y =r A -M m -Y»+Ar r M-Y t -(SM-^Y (8) 

Whenever S( = ( c - Cm = 0, then SM = M c - M m = 
0, Ar^ = 0. Let again Ao be the smaller eigenvalue (in 
absolute value). If Ao ^ it is possible to check that 
([5]) can be fulfilled by setting S( = 0: for regular points, 
transport codes can potentially attain a perfect recon- 
struction of measured value, provided the set of trial pro- 
files includes the solution given by the matricial method. 
This is no longer true when Ao = 0, since in this case 
the l.h.s. of © is aligned along the other eigenvector Ei, 
whereas the r.h.s. has generically components along both 
directions. At a singular point a transport code provides 
a natural regularization because of the limited flexibility 
of the set of Y profiles available for the minimization. In 
terms of Eq. ||HJ) the code generates the minimum finite 
value of 5((r s ) over the set of trial profiles providing a 
finite value of <5M^(r s ) which cancels the Eo component 
of Sr(r s ), that exists for 6((r s ) = 0. More precisely the 
choice of all the <5C(r s )'s is done simultaneously, because 
of the global norm used by the codes. Therefore this 
optimized finite value of SC{r s ) may not be the smallest 
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one. It is independent of <5*Y(r s ) if (5M<5£(r s ) is small. 
Therefore S((r a ) does not measure the accuracy of the 
calculation of Y(r s ). Since 5M(5C(r s ) has a random sign, 
it does not improve this accuracy, even if it is not negli- 
gible. Finally the finite value of S((r s ) also modifies the 
component of Y(r s ) along Ei, which increases the un- 
certainty of the calculation. 

Summarizing, the matricial approach: (I) highlights that 
the source term plays a critical role for profile calculation; 
the best situation is obtained when the source is local- 
ized at the plasma outer edge: then, the inverse problem 
is well-conditioned and the error on the source estimate 
does not enter the calculation of the profile for smaller 
radii. The worst case is conversely expected to be when 
the source is spred radially. Then the calculated profile 
may be strongly influenced by the somewhat arbitrary 
regularization. (II) Points out that the common strat- 
egy of least squares minimization performed by transport 
codes actually hides some subtle practicalities. 
So far, we have implicitly assumed that measurements 
are performed on a very fine mesh of points, such that £ 
may be treated as a continuous variable. In real cases, 
measurements are taken on a discrete and often quite 



coarse mesh, 



1,,..,N. This raises issues of 



interpolation and extrapolation, needed to compute the 
derivatives and the integrals involved in Eq. ((4]), and em- 
phasizes the intrinsic reverse character of our approach 
with respect to transport codes solving Eq. ((T|): the ma- 
tricial approach computes Y only at the discrete set of 
measurement points rj , but needs a continuous interpola- 
tion of the data A and <f> in order to compute the deriva- 
tives involved in Eq. Conversely, transport codes do 
need continuous profiles for {x, V) and only the knowl- 
edge of A, 4> at discrete points. 

From Eq. ([5]) one can estimate the uncertainty associ- 
ated to the calculation of Y(r). Y depends on 3 x N 
experimental quantities {A{\ , {<fii} , {Si}. We label col- 
lectively these quantities as 1?^, k = 1, 37V. Let S'&k be 
an estimate of the associated uncertainty and dk = gf- ■ 
Thus: 



SY = E 

k.l 

+ E, 



r-E, 

A, 



A; 



r-d fc E, 

Ai 



r-E, 



(9) 



Generally, the Sdk 's should be treated as stochastic vari- 
ables picked up from some statistical distribution. Ac- 
cordingly the vector 6Y(rj) spans an area W around the 
point Y(fj): the uncertainty on this quantity is geomet- 
rically displayed as a roughly elliptic region with low ec- 
centricity if Ao/Ai w O(l) . Near the points where Ao is 
small the term proportional to Aq 2 dominates: 



5Y 



-Ei 



r-E ^ 



\2 



A; 



d 1 ? fc (9 fc Ao)+E 1 



Thus W is stretched along Eo, as qualitatively assessed 
earlier. 



Till now we considered a single experimental setup en- 
dowed with a single modulation frequency and source. 
However it is easy to generalize to experiments wherc- 
under the same transport conditions-multiple harmon- 
ics are measured, or the location of sources is changed. 
Then, the above procedure may be repeated, and the in- 
dependent domains of uncertainty W compared. If for 
all points they have a non vanishing intersection, this 
should provide a smaller global uncertainty, and thus im- 
prove the accuracy of the profile calculation. If one of 
the intersections vanishes, the assumptions of the calcu- 
lation must be questioned. In particular one may wonder 
whether: (I) the uncertainties have been underestimated 
and whether the correct W's are larger than estimated; 
(II) Some source terms have not been computed correctly 
or are missing; (III) It is not true that the independent 
measurements correspond to the same (x(r),V(r)) pro- 
file; (IV) Eq. ([TJ fails: transport is not of the advection- 
diffusion type. Whenever the modulated source has sev- 
eral harmonic components, by virtue of the linearity of 
Eq. ([2]), each harmonic can be treated as a separate mea- 
surement. This may improve the accuracy in the regions 
where Ao is not small. If Ao is small, since Eo has almost 
the same orientation for all harmonics, the precision of 
the calculation cannot be improved. 
In order to make visual the above statements we pro- 
pose here below a check using synthetic data. Spatial 
profiles of transport coefficients and sources are given in 
advance (Fig. la,b) and used in Eq. ([2]) to produce syn- 
thetic (A, </>) datasets (Fig. lc,d). Notice that the central 
source S2 produces a singular point close to half radius 
(i.e., near the location of its peak), whereas the edge 
source S± does not. The matricial method is then em- 
ployed on these data to check that the original transport 
coefficients are correctly recovered back (Fig. le,f). Fi- 
nally, we add a finite uncertainty to the " measurements" 
in order to mimic experimental error: to each point is 
added a random contribution taken from a normal dis- 
tribution with a mean amplitude equal to 1% (both in 
the amplitude and the phase). The new datasets have 
been then smoothed using a moving average, in order 
to avoid too brusque variations, that would deteriorate 
the quality of the derivatives, and were interpolated using 
3rd-order splines. Then the coefficients (x, V) are recom- 
puted, repeating the whole procedure for a total of 400 
independent statistical runs at x = 0.41, i.e., close to the 
singularity for the source S2. Figs. (lg,h) show how the 
different estimates spread around the true value. For the 
" regular" case Si all the estimates have a relatively small 
spread, unlike S2, whose data align along the eigenvector 
Eo, as predicted earlier, when Ao — > 0. The width of the 
spread is remarkable, which stresses again the ill-posed 
nature of the problem. 

To summarize, the matricial approach (MA) is very 
light computationally. Indeed it avoids the heavy spatio- 
temporal integration of Eqns. (|H2p and the iterative 
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FIG. 1: (a) Profiles of \ and V; (b) Profiles of the two sources, 
Si , S2 • Both sources are modulated with frequency v = 1. 
(c) Amplitude A of the signal produced by solving Eq. ((2| . 
Circles refer to source Si, diamonds to 52- (d) Phase <f> of 
the same signal, (e) Symbols are the transport coefficients 
computed backwards from Eq. Q, superimposed to the true 
transport coefficients, for the simulation with source Si. (f) 
The same for source 5*2. (g) Scatter of V) from Monte- 
Carlo simulation using source Si at x = 0.41. The black 
segment is parallel to the local eigenvector Eq. (h) The same 
as figure (g) for data from source 5*2. 



least-squares minimization procedure over the set of trial 
functions. The MA provides a clear geometrical foun- 
dation to the nature and size of uncertainties in profile 
reconstruction. The reconstruction radius-by-radius en- 
ables to see how different are the uncertainties over Y 
as a function of r and to correlate them with the pres- 
ence of a source. The MA yields a higher precision in 
the reconstruction of transport profiles than transport 
codes, provided the uncertainty on the estimate of the 
derivatives of A and <fi is not high. Indeed the MA is not 
restricted by the a-priori guess of the trial profiles, but 
by that of these derivatives, which is much more reliable 
and controllable. It may provide an assessment of profile 
predictions already done with transport codes. A poste- 
riori some regularization may be applied to MA results. 
When Ao is small at a given radius r s , the regularization 
needed to provide a reasonable estimate of Y(r s ) can be 
defined in an explicit way, while this is only implicit in 
the transport code approach. For instance, if a singular 



point r s is in between two nearby regular ones r, and fi+i, 
one may require Y(r s ) to be on the straight line joining 
Y(fj) and Y^j+i). Overlapping the uncertainty inter- 
vals for various experiments where the same transport is 
assumed to hold, provides either a way to improve the 
precision of the reconstruction (case of a non vanishing 
overlap) or to prove the set of initial assumptions in the 
reconstruction to be wrong (case of a vanishing overlap) . 
The MA can help designing a priori the combination of 
perturbation measurements susceptible of improving the 
precision of the reconstruction of transport profiles. The 
MA requires a single boundary condition only: the van- 
ishing flux at r = 0, which is a rigorous constraint. In 
contrast the calculation via Fokker-Planck equation of 
the CiW s requires a second outer boundary condition 
which is generally known with some uncertainty. Very 
often this condition is provided by the data of the outer- 
most chord. Then the matricial approach has the benefit 
to keep the outermost chord data available for the pro- 
file calculation. It also avoids the uncertainty included in 
this data to contaminate the profile calculation at smaller 
radii. 
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